function [ffactor, A0PF, Ltilde] = penalty_fun_simple(Sigma,n,p,F,J,Hor,nperPF)
% function [ffactor, A0PF, Ltilde] = penalty_fun_simple(Sigma,n,p,F,J,Hor,nperPF)
%-----------------------------------------------------------------------------
% Compute Impulse Responses and A0 matrices for
% the penalty function approach under a simple fiscal rule.
%-----------------------------------------------------------------------------
% Inputs:
%   Sigma: (draw of) covariance matrix
%   n: number of variables
%   p: number of lags
%   F: companion matrix
%   J: shock linear combination selector matrix
%   Hor: IRF horizonfor plot
%   nperPF: period restrictions
%
% Outputs:
%   ffactor: [nvar by n_shocks] impact matrix for identified shocks
%   A0PF:    [nvar by nvar] A0 matrix from sign restrictions
%   Ltilde:  [Hor by nvar by n_shocks] IRF matrix
%
% Based on replication codes for Dario Caldara and Christophe Kamps (2017), The Analytics of SVARs: A
% Unified Framework to Measure Fiscal Multipliers, Review of Economic Studies (2017) 84, 1015–1040
%
% Copyright: 2017 Dario Caldara and Christophe Kamps
% Copyright: 2020-2023 Benjamin Born, Francesco D'Ascanio, Gernot J. Mueller, Johannes Pfeifer

%--------------
% Preliminaries
%--------------

global IR;
global ssigma;
global nper;
LC =chol(Sigma,'lower');
A0 = (LC')\eye(size(LC,1));

Omega1 = [LC;zeros((p-1)*n,size(LC,2))];
LtildeTemp    = vm_irf(F,J,LC,4,n,Omega1); %do IRF on lower Cholesky factor to get possible combinations
s = sqrt(diag(Sigma));
ssigma = s;
nper = nperPF;
IR = LtildeTemp;

%Identify forecast shock
temp=Sigma;
temp=temp(:,6:-1:1);
temp=temp(6:-1:1,:); %how would impact look if FE were ordered first
UC_permuted=chol(temp,'lower');

Q_FE=LC\UC_permuted(6:-1:1,1); %compute associated rotation if instead ordered last
%     check=LC*Q_FE;
%     check(end).^2-Sigma(end,end)
%------------------------------
% Identify business cycle shock
%------------------------------

Aeq=Q_FE';
beq=0;
q1ga=rand(n,1); %start for rotation vector
[Qhat1,~] = fmincon(@penalty_bc_shock,q1ga,[],[],Aeq,beq,[],[],@mycon,... %mycon restricts length of rotation vector to 1
    optimset('MaxFunEvals',40000,'MaxIter',20000,'Display','off','Algorithm','active-set'));

%---------------------------------------------------------------------------------
% Identify spending shock (imposing a simple fiscal rule through zero restrictions)
%---------------------------------------------------------------------------------

Z = zeros(n-2,n);
Z(1,1) = 1;
Z(2:3,4:5) = eye(2); % Zero restrictions on elasticities in policy rule (6);
Rx = [Z*A0; Qhat1'; Q_FE'];
%Z*A0: only allow for non-zero contemporaneous elasticity for Y and G;
%combines "simple rule in which the fiscal variable can respond
%contemporaneously only to output" (p. 1019 bottom) with "restrict the
%contemporaneous interaction between the two fiscal instruments to
%zero" (p. 1024)
%Qhat1: imposes orthogonality to BC shock, see (6) on p. 91 of M/U (2009)
Nx=null(Rx); %get basis for null space, i.e. set of solution to Rx=0
Xhat=randn(n,1);
Qhat2 =Nx*Nx'*Xhat/norm(Nx'*Xhat);

%---------------------------------------------------------------------------------
% Identify Tax shock  (imposing a simple fiscal rule through zero restrictions)
%---------------------------------------------------------------------------------

Z = zeros(n-2,n); Z(1,2) = 1; Z(2:3,4:5) = eye(2); % Zero restrictions;
Rx = [Z*A0; Qhat1'; Q_FE'];
%Z*A0: only allow for non-zero contemporaneous elasticity for Y and T
%Qhat1: imposes orthogonality to BC shock, see (6) on p. 91 of M/U (2009)
Nx=null(Rx);
Xhat=randn(n,1);
Qhat3 =Nx*Nx'*Xhat/norm(Nx'*Xhat);

%---------------------------------------------
% Calculate elasticities and impulse responses
%---------------------------------------------

Q = [Qhat1 Qhat2 Qhat3 Q_FE];
ffactor = LC*Q;
A0PF    = A0*Q;

Omega1 = [ffactor;zeros((p-1)*n,size(ffactor,2))];
Ltilde  = vm_irf(F,J,ffactor,Hor,n,Omega1);

end